#ifndef __PARSER_H__
#define __PARSER_H__

#include <jsoncpp/json/json.h>
#include <filesystem>
#include <iostream>
#include <vector>
#include <fstream>
#include <filesystem>
#include "MultiGrid.h"
#include "Functions.h"
#include "Timer.h"

#define SECTOR "-----------------------------------"

using namespace std;
namespace fs = std::filesystem;

const string input_dir = "./inputs/";

GridOp::Linear_Interp<Real, 1> Linear_1d;
GridOp::Linear_Interp<Real, 2> Linear_2d;
GridOp::Quadra_Interp<Real, 1> Quadra_1d;
GridOp::Quadra_Interp<Real, 2> Quadra_2d;
GridOp::Injection<Real, 1> Injec_1d;
GridOp::Injection<Real, 2> Injec_2d;
GridOp::FullWeight<Real, 1> FullWeight_1d;
GridOp::FullWeight<Real, 2> FullWeight_2d;

inline void prtNorm(std::ostream &os, double *v)
{
    for (int i = 0; i < 4; ++i)
        os << v[i] << ' ';
    os << endl;
}

MultiGridSolver<1> parse1dSolver(Json::Value &root)
{
    printf("Solving 1-dimensional problem \n");

    GridOp::Interpolation<Real, 1> *Interp_1d;
    std::string interp_type = root["Interpolation"].asString();
    std::cout << "Interpolation Operator : " << interp_type << std::endl;
    if (interp_type == "linear")
        Interp_1d = &Linear_1d;
    else if (interp_type == "quadratic")
        Interp_1d = &Quadra_1d;

    GridOp::Restriction<Real, 1> *Restric_1d;
    std::string restric_type = root["Restriction"].asString();
    std::cout << "Restrcition Operator : " << restric_type << std::endl;
    if (restric_type == "injection")
        Restric_1d = &Injec_1d;
    else if (restric_type == "full-weighting")
        Restric_1d = &FullWeight_1d;

    MultiGridSolver<1> solver(Interp_1d, Restric_1d);

    int set_idx = root["FunctionSet"].asInt();
    std::string type = root["BoundaryCondition"]["type"].asString();
    std::cout << "Solving BVP with boundary condition : " << type << std::endl;
    return std::move(solver);
}

MultiGridSolver<2> parse2dSolver(Json::Value &root)
{
    printf("Solving 2-dimensional problem \n");

    GridOp::Interpolation<Real, 2> *Interp_2d;
    std::string interp_type = root["Interpolation"].asString();
    std::cout << "Interpolation Operator : " << interp_type << std::endl;
    if (interp_type == "linear")
        Interp_2d = &Linear_2d;
    else if (interp_type == "quadratic")
        Interp_2d = &Quadra_2d;

    GridOp::Restriction<Real, 2> *Restric_2d;
    std::string restric_type = root["Restriction"].asString();
    std::cout << "Restrcition Operator : " << restric_type << std::endl;
    if (restric_type == "injection")
        Restric_2d = &Injec_2d;
    else if (restric_type == "full-weighting")
        Restric_2d = &FullWeight_2d;

    MultiGridSolver<2> solver(Interp_2d, Restric_2d);

    int set_idx = root["FunctionSet"].asInt();
    std::string type = root["BoundaryCondition"]["type"].asString();
    std::cout << "Solving BVP with boundary condition : " << type << std::endl;
    return std::move(solver);
}

void SolveBVP2d(const string &filename, const string &relative = "")
{

    std::cout << SECTOR << std::endl;
    std::cout << "Solving the problem specified by : " << filename << std::endl;

    string prefix = filename.substr(0, filename.size() - 5);
    const string out_dir(relative + "./out/2D/" + prefix + "/");
    fs::create_directories(relative + "./out/2D/" + prefix);

    ifstream test_file(relative + input_dir + "2D/" + filename);
    if (test_file)
    {
        Json::Value root;
        Json::Reader reader;
        if (!reader.parse(test_file, root))
        {
            std::cout << "Invalid input stream!" << std::endl;
            return;
        }

        int set_idx = root["FunctionSet"].asInt();
        Function<2> u;
        u = funcs_2d[set_idx]->u;
        double a = root["BoundaryCondition"]["alpha"].asDouble(),
               b = root["BoundaryCondition"]["beta"].asDouble();
        MultiGridSolver<2> solver = parse2dSolver(root);
        FunctionSet<2> *funcs = funcs_2d[set_idx];
        BoundaryConditions<2> BC(*funcs, a, b);
        Function<2> g[4];
        for (int i = 0; i < 4; ++i)
        {
            BC.bounds[i].g = [=](Point<2> p)
            { return a * u(p) + b * (funcs->g_partial[i](p)); };
        }
        solver.define_boundary(std::move(BC));
        solver.assignF(funcs->f);

        string method = root["Cycle"].asString();
        std::cout << "Using Cycle : " << method << std::endl;

        double norm1[4], norm2[4], norm_inf[4];

        // vector<Real> solution;
        for (int i = 0, grid_size = 16; grid_size <= 128; ++i, grid_size *= 2)
        {
            ofstream data(out_dir + to_string(grid_size) + ".txt");
            solver.discretization(grid_size);
            vector<Real> solution(pow(grid_size + 1, 2));
            vector<Real> LU_solution(pow(grid_size + 1, 2));

            if (grid_size <= 64)
            {
                Timer LU("LU");
                solver.solve(solution, "LU");
            }

            {
                Timer MG(method);
                solver.solve(solution, method, 2, 2);
            }

            for (auto p : solver.getGrid())
            {
                data << p[0] << ' ';
            }
            data << std::endl;
            for (auto p : solver.getGrid())
            {
                data << p[1] << ' ';
            }
            data << std::endl;

            vector<double> r;

            int grid_num = solver.getGrid().size();

            double c = solution[0] - u(solver.getGrid()[0]); // for neumann condition

            for (int j = 0; j < solution.size(); ++j)
            {
                r.push_back(solution[j] - u(solver.getGrid()[j]) - c);

                data << solution[j] - c << ' ';

                norm1[i] += fabs(r[j]);
                norm2[i] += r[j] * r[j];
                norm_inf[i] = std::max(norm_inf[i], fabs(r[j]));
            }
            norm1[i] /= grid_num;
            norm2[i] = sqrt(norm2[i] / grid_num);

            data << std::endl;
            std::cout << "n : " << grid_size << std::endl;
            std::cout << "1-norm : " << norm1[i] << std::endl;
            std::cout << "2-norm : " << norm2[i] << std::endl;
            std::cout << "infinity-norm : " << norm_inf[i] << std::endl;
            std::cout << std::endl;
        }
        ofstream norm(out_dir + "norm.txt");

        prtNorm(norm, norm1);
        prtNorm(norm, norm2);
        prtNorm(norm, norm_inf);

        std::cout << "Convergence rate: " << std::endl;
        std::cout << "1-norm : " << my_conv_rate(norm1, 4) << std::endl;
        std::cout << "2-norm : " << my_conv_rate(norm2, 4) << std::endl;
        std::cout << "infinity-norm : " << my_conv_rate(norm_inf, 4) << std::endl;
    }
    else
        std::cout << "Failed to open the json file!" << std::endl;
}

void SolveBVP1d(const string &filename, const string &relative = "")
{

    std::cout << SECTOR << std::endl;
    std::cout << "Solving the problem specified by : " << filename << std::endl;

    string prefix = filename.substr(0, filename.size() - 5);
    const string out_dir(relative + "./out/1D/" + prefix + "/");
    fs::create_directories(relative + "./out/1D/" + prefix);

    ifstream test_file(relative + input_dir + "1D/" + filename);
    if (test_file)
    {
        Json::Value root;
        Json::Reader reader;
        if (!reader.parse(test_file, root))
        {
            std::cout << "Invalid input stream!" << std::endl;
            return;
        }

        int set_idx = root["FunctionSet"].asInt();
        Function<1> u;
        u = funcs_1d[set_idx]->u;
        double a = root["BoundaryCondition"]["alpha"].asDouble(),
               b = root["BoundaryCondition"]["beta"].asDouble();

        MultiGridSolver<1> solver = parse1dSolver(root);
        FunctionSet<1> *funcs = funcs_1d[set_idx];
        BoundaryConditions<1> BC(*funcs, a, b);
        for (int i = 0; i < 2; ++i)
        {
            BC.bounds[i].g = [=](Point<1> p)
            { return a * u(p) + b * (funcs->g_partial[i](p)); };
        }
        solver.define_boundary(std::move(BC));
        solver.assignF(funcs->f);

        string method = root["Cycle"].asString();
        std::cout << "Using Cycle : " << method << std::endl;

        double norm1[4], norm2[4], norm_inf[4];

        vector<double> solution;
        for (int i = 0, grid_size = 32; grid_size <= 256; ++i, grid_size *= 2)
        {
            solver.discretization(grid_size);
            ofstream data(out_dir + to_string(grid_size) + ".txt");
            solution.resize(grid_size + 1);
            solver.solve(solution, method, 2, 2);

            for (auto p : solver.getGrid())
            {
                data << p[0] << ' ';
            }
            data << std::endl;

            vector<double> r;

            int grid_num = solver.getGrid().size();

            double c = solution[0] - u(solver.getGrid()[0]); // for neumann condition

            for (int j = 0; j < solution.size(); ++j)
            {
                r.push_back(solution[j] - u(solver.getGrid()[j]) - c);

                data << solution[j] - c << ' ';

                norm1[i] += fabs(r[j]);
                norm2[i] += r[j] * r[j];
                norm_inf[i] = std::max(norm_inf[i], fabs(r[j]));
            }
            norm1[i] /= grid_num;
            norm2[i] = sqrt(norm2[i] / grid_num);

            data << std::endl;
            std::cout << "n : " << grid_size << std::endl;
            std::cout << "1-norm : " << norm1[i] << std::endl;
            std::cout << "2-norm : " << norm2[i] << std::endl;
            std::cout << "infinity-norm : " << norm_inf[i] << std::endl;
            std::cout << std::endl;
        }
        ofstream norm(out_dir + "norm.txt");

        prtNorm(norm, norm1);
        prtNorm(norm, norm2);
        prtNorm(norm, norm_inf);

        std::cout << "Convergence rate: " << std::endl;
        std::cout << "1-norm : " << my_conv_rate(norm1, 4) << std::endl;
        std::cout << "2-norm : " << my_conv_rate(norm2, 4) << std::endl;
        std::cout << "infinity-norm : " << my_conv_rate(norm_inf, 4) << std::endl;
    }
    else
        std::cout << "Failed to open the json file!" << std::endl;
}

void epsil_test(const string &relative = "")
{
    string filename = "FQV.json";
    std::cout << SECTOR << std::endl;
    std::cout << "Solving the problem specified by : " << filename << std::endl;

    string prefix = filename.substr(0, filename.size() - 5);
    const string out_dir(relative + "./out/1D/" + prefix + "/");
    fs::create_directories(relative + "./out/1D/" + prefix);

    ifstream test_file(relative + input_dir + "1D/" + filename);
    if (test_file)
    {
        Json::Value root;
        Json::Reader reader;
        if (!reader.parse(test_file, root))
        {
            std::cout << "Invalid input stream!" << std::endl;
            return;
        }

        int set_idx = root["FunctionSet"].asInt();
        Function<1> u;
        u = funcs_1d[set_idx]->u;
        double a = root["BoundaryCondition"]["alpha"].asDouble(),
               b = root["BoundaryCondition"]["beta"].asDouble();

        MultiGridSolver<1> solver = parse1dSolver(root);
        FunctionSet<1> *funcs = funcs_1d[set_idx];
        BoundaryConditions<1> BC(*funcs, a, b);
        for (int i = 0; i < 2; ++i)
        {
            BC.bounds[i].g = [=](Point<1> p)
            { return a * u(p) + b * (funcs->g_partial[i](p)); };
        }
        solver.define_boundary(std::move(BC));
        solver.assignF(funcs->f);

        string method = root["Cycle"].asString();
        std::cout << "Using Cycle : " << method << std::endl;

        // double norm1[4], norm2[4], norm_inf[4];
        std::cout << "Begin epsilon test ....." << std::endl;

        for (int n = 32; n <= 256; n *= 2)
        {
            std::cout << "Using grid size : " << n << std::endl;
            solver.discretization(n);
            vector<double> solution(n + 1);
            double epsil;
            for (epsil = 1e-8; epsil >= 2.2e-16; epsil /= 2)
            {
                solver.solve(solution, "VCycle", 2, 2, epsil, 100);
                if (solver.getResidNorm() < epsil)
                    std::cout << "Epsilon : " << epsil << " has been reached." << std::endl;
                else
                {
                    std::cout << "Failed to reach epsilon : " << epsil << std::endl;
                    break;
                }
            }
            std::cout << "The minimum reachable epsil : " << epsil << std::endl;
            std::cout << std::endl;
        }
    }
    else
        std::cout << "Failed to open the json file!" << std::endl;
}

#endif // __PARSER_H__